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ABSTRACT 

The properties of dense granular systems are analyzed from a hydrodynamical point of view, based on conservation 
laws for the particle number density and linear momentum. We discuss averaging problems associated with the nature 
of such systems and the peculiarities of the sources of noise. We perform a quantitative study by combining analytical 
methods and numerical results obtained by ensemble-averaging of data on creep during compaction and molecular 
dynamics simulations of convective flow. We show that numerical integration of the hydrodynamic equations gives 
the expected evolution for the time-dependent fields. 
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,0, : 1. INTRODUCTION 

The interest in granular materials goes back to the early sixties, when industries needing to process powders required 
' ^ ■ more and more control over the quality of their products. One may imagine a variety of shapes for mixers and many 
I ■ different ways to drive the grains inside into motion; a down-to-earth approach would consist of an answer for the 
' question: in which of these mixers should I spend money? This can be done if one can compare their performances 
, by means of realistic models that reproduce true flow patterns. 

■ 

, In order to answer aspects of this problem, and based on the gas-like appearance of the compounding particles, 

fS| ' attempts were made to apply a hydrodynamical description to granular systems. It has been more than a deca " 
, since the pioneering contributions to the understanding of granular materials from a hydrodynamical point of viewl 
This understanding, however, is still far from being complete. We will address here several questions of theoretical 
and practical interest, namely fluctuations and averaging problems intrinsic to granular nature, modelling of dense 
Ch ' granular flows by means of hydrodynamic equations and boundary conditions, comparisons with experimental data 
I and results from molecular dynamics simulations, and fluctuations and mixing problems. 

c : 

O . 2. THE SOURCES OF FLUCTUATIONS 

^ The behavior of a granular mass is fundamentally different from that of typical fluids. For instance, consider ordinary 
hydrodynamic fluctuations: the lack of an intermediate length scale -much larger than the typical diameter of the 
grains but much smaller than the size of the system, makes the thermodynamic limit unnatainable and therefore 
hydrodynamic fluctuations may subsist in the continuum limit. In addition, the fact that solid grains in a dense 
arrangement can't be regarded as points in any length scale leads to a second source of fluctuations; it operates 
at distances much larger than the typical diameter of the particles and is related to the appearance of extensive 
arrangements inside the system. Given a mean density not far away from the close-packing limit, there is a fraction 
of the free volume which is still available and can be distributed in many ways. Some of these configurations will flow, 
some will not allow net motion, so different actual realizations in the configuration space of the particles may lead to 
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departing dynamic properties and, ultimately, in generating completely different time-sequences. This contribution 

is what we call non-local noise. Both, local and non-local noise, are always present but their relative importance 
strongly depends on the forcing applied to the system, measured by the parameter J^g = f / pg, f being the volume 
density of the forcing and g the acceleration of gravity. 

Therefore, one can distinguish two sources of fluctuations and two types of averaging, one over local noise and 
the other over different configurations or realizations. 

We propose the existence of two regimes; in the weak- forcing limit, J-'g ^ 1, the non-local component of the 
noise dominates, and consequently one expects very long relaxation rates and non-self-averaging quantities. In this 
limit, only ensemble averaging is meaningful, as in every possible realization the system explores a small region of 
the configurational space. The glass-like behavior is most apparent. In the opposite limit, !Fg ^ 1, the system is 
not easily trapped in an immobile arrangement, and one can safely assume that in a sufficiently long time it explores 
the representative part of the configuration space. Time averaging can substitute ensemble averaging only in this 
limit. Consistently with this picture, the critical density, pc, is not unique in general, but is a distributed quantity 
depending on the configurational state, F. Actually, we shall see that experimental data on compaction at weak 
forcing display quenched behavior, and the final density may vary over more than 10%. Instead, our numerical study 
of granular convection under strong forcing indicates a much narrower histogram of maximal achieved densities, Pc, 
despite the fact that the number of particles and the used number of samples for averaging were much smaller. 
Consider that for high forcing, the mean density evolves within a comparatively wide margin and the system very 
rarely, if ever, visits any of the corresponding trapping configurations. In the mean-field limit, one can consider pc 
as an unique constant, which would correspond to the strict close-packing limit. 

In spite of the complexity of this picture, we show that it is possible to understand both weak and strong forcing 
limits within the frame of hydrodynamic equations, which are in general stochastic equations including both local 
and non-local noise. It is beyond of the scope of this paper a detailed analysis of the properties of such equations, 
which is extensively done elsewhere. Our aim is to present some results of our study of the evolution of the mean 
hydrodynamic fields, comparing them with their observed behavior in a sample of cases. 



3. HYDRODYNAMIC EQUATIONS 

By hydrodynamic equations we mean balance equations for mean mesoscopic hydrodynamic fields. We will focus on 
the equations for the conservation of mass and linear momentum, leaving the energy balance for a later discussion. 
As said in the previous section, the lack of an intermediate length scale -contrary to what happens in simple fluids 
for instance, containing a sufficiently large number of particles such that local fluctuations fade away, adds on the 
above exposed problem of non-local fluctuations due to the intrinsic granular nature of the system, operating at 
length scales where hydrodynamic fields can already be defined. The former can be modelled as the usual additive 
stochastic contribution to the dissipative flows and comes related to the existence of some kind of "temperature", 
the latter enters via distributed kinetic coefficients (depending on the conflgurational state F). In the continuum 
approach, conservation of mass and linear momentum read, in the Stokes approximation, 

dtP + W-pv = 0, (1) 

Pdm = ^ + U, (2) 

where /» are the components of the volume density of forcing. Note that these equations express very basic laws. 
However, they are valid in the usual form when the averages decouple, that is, when the average of the "microscopic" 
linear momentum equals the product of the averages of p and v. However, due to the mesoscopic natiire of the 
averaging procedure, this is not guaranteed. It will work if the local fluctuations are small enough. We shall discuss 
the fact that local fluctuations of the velocity fleld are very small in dense granular flows, except the very few moments 
when the granular mass undergoes dramatic changes, and close to the boundaries of the system. 

As for the terms constituting the stress tensor cry a few comments are in order. Provided that the granular 
particles may be modelled as sufficiently hard spheres, we neglect any elastic contribution other than that introduced 



by pressure effects, and assume tfiat the non-equilibrium part of is a local functional of the derivatives of the 
local velocity, 

(3v ' 3v 2 \ 
^ + ^ - 3% V • vj + C(/9, r)% V • V + . (3) 

where p is the pressure and 77 and C the shear and volume viscosities, respectively. 

3.1. The role of temperature and the pressure term 

Our molecular dynamic simulations indicate that the evolution of the velocity field in dense granular fiow is nearly 
independent from the granular temperature, defined as the mean fluctuational part of the velocity. Effectively, well 
inside the bulk, the quantity Sv/v, measuring such fiuctuational deviations from the mean velocity v, is typically 
about six orders of magnitude smaller than close to the boundary. This and other evidences allow us to suppose that 
granular flows in dense systems can't be sustained by a temperature-based mechanism alone -unlike Rayleigh-Benard 
convection in simple fluids, for instance. Note that the equation for the energy balance has been omitted; consistently 
with the observation that the granular temperature plays no significant role, its evolution appears decoupled from 
the previous system of equations^. Similarly, one cannot account for elastic contributions in the high density limit 
only by using a thermal pressure. Therefore, one has to model such terms by means of an artificial equation of state 
which must help to resolve the delicate limit p Pc- We assumed the most simple dependence, p — Po/(l ~ p/Pc), 
where pq represents a certain constant, in our numerical integration of the hydrodynamic equations in the strong 
forcing regime. 

3.2. The viscosity 

Accordingly, the model that we adopt for the viscosity is not thermal, but glass-like. In dense clusters, in order 
to move, a complex rearrangement of particles has to occur making use of voids. Similar properties are exhibited by 
glasses. Available experimental data and our numerical results indicate the presence of a factor exp[c/(l — p/pc))] 
in the jpean flow rates, where c is a dimensionless number. This formula is related to the Vogel-Fultcher law for 
glasseaa; it measures the number of attempts needed for one step in the direction of average flow in a dense granular 
system. Sufficiently close to pc, we then expect a shear viscosity of the form 

,(p,r)=,o(p,r)exp(^— ^). (4) 
and a similar dependence for the bulk viscosity, F). 

4. TEST OF THE HYDRODYNAMIC EQUATIONS 

We discuss some representative examples of each regime. For the weak forcing limit, available experimental data 
on compaction of sand during tapping experimentsQ provide the necessary reference. For the strong forcing limit 
we perform extensive ensemble averaging of samples generated by molecular dynamic (MD) simulations of vertical 
and horizontal shaking, comparing the resulting hydrodynamic fields with those generated integrating numerically 
the hydrodynamic equations. We shall also show analytical results on the cycle-averaged velocity profiles in vertical 
shaking that fit experimental profiles.El 

* Actually, we observed that the rate of energy release, coming from inellastic collisions, was clearly correlated with the viscous heating 
(see next subsection for the model for the viscosity). These terms are therefore responsible for the evolution of the temperature field, 
since the energy loss is proportional to some power of the granular temperature, T. We checked the functional dependency, resulting in 
the expected T^/^. 



4.1. Weak forcing limit. Application to compaction experiments. 

These results provide additional support for the existence of non-local noise, and evidence that the mean flows are of 
hydrodynamic nature even in very dense limits. Beginning with a loosely packed sand at volume fraction po = 0.57 
the authors report a logarithmic density growth, 

P^~P^'^= 1 + BHl + t/ry 

where A,B,T,pc are four fitting parameters. It can be shown that it is possible to retrieve such a dependence by 
integration of hydrodynamic equations. Omitting further details about calculations and average over non-local noise, 
one finds after integration of the 1-dimensional version of (|l|,||) at late times 

H{J^dtF)/cvo]' 

Pc and c, already averaged over non-local noise, are constants which depend on, say, the amplitude of forcing, but 
do not change over time. The quantity F is related to the integral of the density of forcing and is left unspecified. 
Equation can be compared with experimental fit, (see Fig. ^) assuming Jl^dtF ^t{F), where the average is taken 
over a period of repeated tapping. We find A/B = c, t = cfjo/ (F). This is a three parameter fit. It demonstrates that 
hydrodynamics may be used for analyzing experimental data. Different fitting values oi Pc, c support the assumption 
that granular configurations with different pc, c (different states F) do not communicate at weak forcing. Now wc 
can be more specific in what we mean by "weak" forcing. Note that the fit is satisfactory at late times. At early 
times that situation is different: the higher is the value of the longer it takes before the fit is any good. If we plot 
the density when the deviation from the fit is a few percent versus the forcing parameter we get Fig. |l|b. It suggests 
that the fits are good when the assumption of quenched values of pc and c is fulfilled. Therefore, one may expect 
that the non-local noise is quenched above a certain line in the {p, J-g) plane. 




Figure 1. a) Fits of the experimental curves of the density dependence on the number of taps, from Ref. ^. 
Tg = 1.4, 1.8, 2.3, 2.7 from the bottom to top. The following parameter values were used to fit the curves in the same 
order: (F)/r/o = (6.9 ■ 10-^ 4.1 • 10^, 3.1 ■ 10\ 1.1 ■ IQ-^), c = (0.92, 0.029, 0.18, 1.35), pc = (0.5985, 0.599, 0.62, 0.67). 
b) Region of quenched non-local noise as a function of the density and the forcing parameter Fg. 

4.2. Strong forcing limit. Results from MD simulations. 

This case is object of a more complete study. Consider for example periodic vertical shaking of sand under gravity. It 
is well known that sand in such conditions develops typical convective rolls, with particles going upwards inside, and 
downwards along the vertical walls. The motion is evidenced, for example, by the bulging colored stripes resulting 



from Magnetic Resonance Imaging experiments,El which represent cycle-averaged displacements. Again, it is possible 
to show that the system of equations (|l|,^ can be integrated to give a good fit of the experimental data. We will 
omit here the detailed derivation, which can be found in Ref. |6[ Imposing the geometry of a tall container of width 
L in which the density slightly decreases with height, assuming Eq. (Q) for the viscosity and the incompressibility 
condition (see step 4, below, for the reasons of such approximation), one finds 



, ^ cosh(fcx) 
1 — kL- 



sinh(fcL) 



(7) 



which is solution of Eq. (g) for the cycle-averaged vertical velocity, in a two dimensional geometry. Here represents 
the scale of variation of the viscosity in the linear approximation, k an arbitrary inverse length scale and vq a 
characteristic amplitude of the velocity, unspecified. See Ref. || for experimental evidence of Eq. (^. 

A major understanding of the motion requires, however, a time-resolved analysis, and we show how this can be 
done via molecular dynamics simulations. We will not discuss details about the simulations here, it sufRces to say 
that we used a polydisperse sample of 2000 soft spheres in a rectangular container, and the chosen values for the 
friction parameters reproduce correctly the experimental MRI images mentioned above. J-'g was about 2 in the study 
of vertical shaking, and 9 in the case of horizontal shaking. We choose the first for presentation, although the results 
can be extended to the horizontal shaking, leaving aside some peculiarities which are not worth to comment here. 
The following steps summarize the procedure we used, which is schematized in the chart of Fig. |[ 



Numerical fit of velocity profiles 
Determination of pc, c 



MD simulations 
microscopic quantities 



Time discretization 

Recording of positions 
every 1/100 cycle 



Ensemble averaging 

Generation of liydrodynamic 
fields: density, velocity, 
granular temperature 



Test of liydrodynamic equations 
Numerical integration 



Study of boundary conditions 
Determination of the constant 7]q 



Figure 2. Scheme of the procedure followed for the test of the hydrodynamic equations. 



1. Time discretization. Each period of shaking was divided in an equal number of frames, where positions of 
particles were recorded. 

2. Spatial discretization and averaging. Using a high resolution grid, the container was divided in cells of the size 
of the order of one particle. Time averaging -which can replace ensemble averaging in this case, was performed 
with data of the corresponding frames of more of 100 periods of shaking. 

3. Mean density, velocity and temperature (mean fluctuational velocity) were obtained and displayed. Visu- 
alization was done by means of IDL movies, showing smooth, well behaved fields. The sequences revealed 
unsuspected details about the motion that the granular mass experiences during a cycle, totally hidden when 
a cycle-average is performed (which leads to the typical convective rolls and has confused the sand community 
for long time). As an example, we reproduce in Fig. 3a the horizontal component of the velocity, Vr^. Observe 



that the motion is complex and unexpected, in the sense that one cannot infer from the sequence of pictures of 
Vx (neither from Vy , not shown) the direction of the global motion. The evolution of the granular temperature 
is shown in Fig. ^. Since the darkest regions represent the lowest values, it is obvious that its importance is 
restricted to the region close to the free surface, and in the bulk, to a limited number of frames. 

4. Test of the hydrodynamic equations. By using the mean hydrodynamic fields obtained in 3., the system of 
equations (|l|,^) was checked. Selected frames provided fitting values for c and pc (see Fig. ||a). rjo was found to 
be about 300 cpoise by comparing histograms of the tangential force and the velocity gradient close to the walls 
(Fig. P^), whereas the observation that the flow was mostly divergence-free allowed us to neglect the effects of 

C- 

5. Study of boundary conditions. We obtained effective boundary conditions for the flow that reproduce to some 
extent the assumptions of microscopic friction during collisions, but wc also found that the motion of sand 
along the vertical walls comes accompanied by dramatic periodic changes in the density and the stress (Fig. ||). 

6. The previous results were used to integrate numerically the system of hydrodynamic equations. In Fig. ^ we 
show comparatively the sequence obtained for Vx- Observe that, apart from a qualitative agreement (including 
the order of magnitude of the velocity fields) there is room for improvement. 
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Figure 3. a) (left) Ensemble-averaged horizontal velocity, Vx of granular material in the course of vertical 
shaking. 10 x 10 frames taken during one cycle of oscillation are positioned from left to right and from top to the 
bottom. The frames of highest contrast indicate the collision of the granular system with the bottom wall, with 
compression of the material propagating upwards, b) (right) Evolution obtained by numerical integration of the 
mean- field hydrodynamical equations (|l|,^ with the boundary conditions obtained in step 5. In both pictures dark 
areas represent negative values and white ones positive values of the velocity. 
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Figure 4. 10 x 10 sequence for the granular temperature (ensemble-averaged fluctuational velocity) in one cycle 
of horizontal shaking. The positioning of frames is the same as in Fig. |. White areas indicate high temperature 
regions. By comparing with Fig. ^, observe that the granular mass moves in the opposite direction of the temperature 
gradient. 

5. FLUCTUATIONS AND MIXING 

Consider the motion of a test granular particle in any averaged velocity field generated by the procedure exposed 
above (different container shapes, various types of shaking). This is of practical interest in order to determine, for 
example, the effectiveness of a mixing device. In the absence of granular temperature, the motion of the test particle 
would follow the flux lines. Starting from a given position, x„, after one period a new location x„+i will be reached, 
which is the result of the integration of elementary displacements (here equating the number of frames per cycle 1. 
Fig. 1^ shows the results for vertical and horizontal shaking in the same way as the Magnetic Resonance Imaging.Q 
In this case computer imaging simulates the horizontal and vertical motion of colored bands initially parallel, which 
become distorted in the course of one period. The cycle-averaged motion is then reducible to iterative maps, leading 
to mixing. Such mixing is incomplete since there may exist regions undergoing periodic motion, which are not mixed 
at all. It happens, for example, in the center of convection rolls accompanying horizontal shaking, potentially making 
the horizontal shaking of thick layers ineffective for mixing particles. 

The second origin of mixing is due to noise, leading to diffussion which superimposes to the iterative map discussed 
above. The distance to diffuse in order to achieve complete mixing is determined by the largest region containing 
limiting cycles and/or fixed points. Non-local noise at a length scale In and time scale i„ leads to nested displacements 
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Figure 5. a) (left) Fit of tlie ensemble-averaged profiles of vertical velocity, Vy(x), in frames 55-70 (numbered in 
the same way as in Fig. In the fitting procedure the values for the constants c=0.15 and pc= 1.01 max(p) were 
found. Dramatic changes at the boundaries are due to density discontinuities accompanying downward sliding of 
sand along the side walls, b) (right) Comparative histograms of the distribution of tangential force density acting at 
a wall (solid line) and the non-diagonal component of the hydrodynamic stress tensor, am = '?(p)(9wr/97i), where n 
and T stand for normal and tangential directions, respectively. 

correlated in space and time and helps to mix particles at distances exceeding Z„ and times exceeding t„. This subject 
has to be postponed until the properties of non-local noise are specified. 

6. CONCLUSIONS 

1. Hydrodynamic equations provide an adequate theoretical frame for the study of dense granular systems. 

2. Non-local or configurational noise adds on the ordinary hydrodynamic noise, leading to non-selfaveraging 
properties in the limit of dense flows. 

3. Temperature-based mechanisms can be practically dismissed in the description of dense granular flows. 

4. Averaging of data from molecular dynamic simulations is an useful tool to reveal the details of the evolution 
of hydrodynamic fields. 

5. Results of numerical integration of the system of hydrodynamic equations show a qualitative agreement with 
the evolution of hydrodynamic fields. 

6. Study of fluctuations and cycle-averaged motion concern practical applications such as mixing. 
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Figure 6. a) (above) Averaged dependence of granular density close to the vertical walls during one period of 
vertical shaking. The decrease in boundary density starts with the take-off of the granular material from the bottom 
of the container and reaches its maximum just before its landing (maximum acceleration downwards). Different 
symbols indicate the value of the off-diagonal component of the stress tensor, (t^„. There are seven symbols used 
in linear proportion to the increasing value of am- {+, *, •, A, □, x). b) (below) The same data plotted in a 
different representation [dvr/dn vsvr)- Plotting symbols now account for the value of granular density, in the same 
ascending order in the density range 0.55 < p < 0.9. 
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